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We study bistability in the electron transport through a ring of N coupled quantum dots with 
two orbitals in each dot. One orbital is localized (called b orbital) and coupling of the b orbitals 
in any two dots is negligible; the other is delocalized in the plane of the ring (called d orbital), 
due to coupling of the d orbitals in the neighboring dots, as described by a tight-binding model. 
The d orbitals thereby form a band with finite width. The b and d orbitals are connected to the 
source and drain electrodes with a voltage bias V, allowing the electron tunnelling. Tunnelling 
current is calculated by using a nonequilibrium Green function method recently developed to treat 
nanostructures with multiple energy levels. We find a bistable effect in the tunnelling current as a 
function of bias V , when the size N > 50; this effect scales with the size N and becomes sizable 
at N ~ 100. The temperature effect on bistability is also discussed. In comparison, mean-field 
treatment tends to overestimate the bistable effect. 
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I. INTRODUCTION 

Bistability (or hysteresis) has always stood as one 
of the most fascinating phenomena in condensed mat- 
ter physics: its original mechanism can be either mag- 
netic as in novel materials like metal-ion cluster— and 
organic radical crystals 2 -, or electric as in resonant tun- 
nelling diodes^. For resonant tunnelling structures, their 
intrinsic bistability behaviour has continued to attract 
attention^. Recently, a novel optical bistability in polari- 
ton diodes was also reported^. Besides being physically 
intriguing in its own right, bistable characteristic can be 
used to fabricate memory devices^. Nanoscale devices ex- 
hibiting such property can therefore have potential appli- 
cation for low-cost, low-power and high-density memory 
storage. 

Recent theoretical studies have explored the possibil- 
ity of making a memory device of a highly degenerate 
molecular quantum doti or of a single-level molecular 
junction 8 -^; both cases rely critically on strong electron- 
phonon interaction (polaron effect) in order to make 
possible effective attractive electron-electron interaction. 
However, two experiments measuring transport through 
a carbon nanotube quantum dol^ with an eightfold de- 
generate state, and through a single spherical PbSe quan- 
tum dofeii with a sixfold degenerate state did not show 
any bistable characteristics. These results indicate that 
electron-phonon interaction may not be strong enough in 
carbon nanotube quantum dots or PbSe quantum dots 
to induce attractive electron-electron interaction. It re- 
mains inconclusive whether a single quantum dot junc- 
tion can exhibit bistability. 

Recent experiments have demonstrated that 
semiconductor— and metallic*^ quantum dots can 
be chemically engineered into geometrically-ordered 
superlattices (nanocrystals). Nanoscale manipulation 



allows to control the lattice constant and quantum dot 
size with a precision limited by atomic roughness^, 
to make a quantum dot array (QDA) evolve from a 
Coulomb blockade regime to a semiconducting regime 
by tuning the interdot coupling^, and to form a band 
structure from a QDAi£. As a result, QDAs provide a 
well-controlled system for study of strong correlation^ 
as well as a promising integrated electronic device. 

Previous mean-field study-i£ on QDA junctions de- 
scribed by the Falicov-Kim mode l 18 ' 19 predicted bista- 
bility behaviour in current-voltage characteristics, which 
makes devices made of QDAs a good candidate for low- 
power high-density memory devices. Each dot in the 
QDA junction contains localized and delocalized states 
with on-site repulsive Coulomb interaction. Nonequi- 
librium mean-field theory^ produced two coupled tran- 
scendental equations that can be solved for the average 
single-particle occupation numbers as a function of volt- 
age bias. Meta-stable states are then made possible by 
the charge transfer between these two kinds of states, 
leading to bistable tunnelling current. However, in ad- 
dition to doubt that could be raised against mean-field 
treatment in strong Coulomb interaction 2 ^, this simple 
theory can not give any answer to at what size of the 
QDA junction bistability starts to emerge and how it 
scales with the system size — an important issue for ex- 
perimentalists who intend to realize bistability in such 
system and possibly optimize it for future manufactur- 
ing. 

In this paper, to study the size effect, we consider a 
ring of N coupled quantum dots simulated by a multi- 
level Anderson model. The model is tackled by using a 
nonequilibrium Green function metho d 21 ' 22 that is able 
to take into account strong correlations in intra-orbital 
and inter-orbital states. This method was originally de- 
veloped to study a nanostructure tunnel junction with 
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multiple energy levels in the Coulomb blockade regime. 

Our results show qualitative agreement with those ob- 
tained from nonequilibrium mean-field theory: Bistabil- 
ity depends critically on the strength of the Coulomb en- 
ergy U, the tight-binding hopping parameter td, and the 
ratio of the left and right tunnelling rates; the Coulomb 
energy U needs to be large enough for bistability to ap- 
pear. However, our method shows that mean-field theory 
tends to overestimate the bistable effect by incorrectly 
incorporating higher-energy spectral density of the delo- 
calized states, even though the bias is too low to allow 
population at high energy. In particular, we show that 
bistability in the tunnelling current occurs when N > 50 
and becomes sizable around N ~ 100 in a QDA system 
with typical physical parameters; the best td value for 
pronounced bistability is found around {7/4. This size 
estimation serves as an important guide for fabrication 
of QDA junctions aimed at memory devices. Finally, we 
show that increasing temperature tends to destroy bista- 
bility, as does the interdot Coulomb interaction between 
localized states. To preserve bistability, it is then impor- 
tant to keep the localized states well separated, and for 
the temperature to be much smaller than the effective 
tunnelling rates through the localized states. 

The paper is organized as follows: SectionlTIlintroduccs 
the multi-level Anderson model and the approximations 
involved. In Sec lIIII we formulate the tunnelling current 
and derive its simple linear relation with the average oc- 
cupation number. In Sec lIVl we discuss the nonequi- 
librium Green function method. In Sec|Vl we present 
numerical results. We conclude in the end. 



II. 



MULTI-LEVEL ANDERSON MODEL 



To observe bistability in the tunnelling current through 
a QDA, an experimental setupii was proposed, see 
Fig. [Tic). The QDA is embedded in an insulated layer 
sandwiched by two metallic leads. To these two leads 
a voltage bias is applied and current is measured be- 
tween them. The system can be either 2D or ID QDA, as 
shown in Fig. QJa) and Fig. [lib) f° r finite-sized 2D and 
ID QDAs, respectively. Possible QDA candidates are 
molecular solids and nanocrystals that potentially can 
display bulk properties (e.g. band formation). 

To study the size effect, we consider a ID ring-like ar- 
ray of TV coupled quantum dots (see Fig. QJb)). Each QD 
is in fact a combination of a small dot (which hosts a lo- 
calized orbital, called b orbital) and a larger dot (which 
hosts a delocalized orbital, called d orbital), e.g. in a 
core-shell like structure. We then consider the case that 
the b orbital is either s-like or p z -like, and the d orbital 
is p x -like or p y -like, so they have different parities in the 
plane of QDA. The d orbitals of two neighboring dots are 
coupled via a tight-binding interaction t n d- The d and b 
orbitals are connected to two metallic leads (labeled by 
a = L, R) with an applied voltage (see Fig. HJc)). Both 
the intradot and interdot Coulomb interactions Hjj and 



H x u are taken into account. We therefore use a gen- 
eralized multi-level Anderson models to describe such 
system: 
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where c^ ko .(c Q ker) is the creation (annihilation) opera- 
tor of a free electron of momentum k and spin a(= ±1) 
in the a- lead with energy e^; &| lcr (6 n(T ) and d\ l(J (d n! y) are 
the creation (annihilation) operators of electrons of spin 
a in the b and d orbitals of the nth quantum dot with en- 
ergies E n bv and E n do , respectively. Their respective num- 
ber operators are n nha = b\ la b ncr and n nda = d^d^. 
U n ji(j,l = d,b) denotes the intradot Coulomb interac- 
tion between two electrons of j and I orbitals in the nth 
dot, while Vji(R nm ) the interdot Coulomb interaction be- 
tween the nth and the mth dots. For the latter, we shall 
adopt a screened Coulomb potential, i.e., Yukawa- type 
potential, which reads 



i j U x ji j~ exp 



Ri. 



where R nm {= |Rn — R-m|) is the distance of the two dots 
and Ri is the nearest-neighbor distance. Lji is the screen- 
ing length. The constant U X ji denotes the Coulomb in- 
teraction between nearest neighboring dots, typically a 
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U x bbGxp[— Ri/Lbb]- Wc then take the tight- binding cou- 
pling to be real and site independent; so t nc i = td- 

In real space, \nda) states are coupled via a tight- 
binding interaction (TjTJ) , which we consider to be strong in 
favor of formation of a band with finite width^. This mo- 
tivated us to first transform the \nda) states into band 
states, which are eigenstates of the rotation operators 
that leave the ring of dots invariant. To this end, we 
define d-state operators in momentum space 
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where 0„ = 2irn/N . We can then rewrite parts of the 
Hamiltonian as 



Ht ^ taiijibabnfjCakcr ^ tadad^ )(T C a p 1 ^( 7 -\- h.C. , (9) 
anker aprjtr 
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FIG. 1: (Color online) Finite-sized (a) 2D and (b) ID QDAs. 
Each dot is a combination of a large dot (green) and a small 
dot (brown), (c) QDA embedded into an insulated matrix 
sandwiched by two metallic leads, (d) Schematic of electrons 
tunnelling through a QDA in the presence of an external bias 
V. The QDA consists of localized b states (discrete black 
lines) and delocalized d states forming a finite band (grey 
area) . 



fraction of Udd or Udb- tnd is the hopping interaction be- 
tween the d orbitals in two neighboring nth and (n+ l)th 
dots, and taenia is the tunnelling matrix element between 
the state |ka) in the a lead and the state \nla) in the nth 
dot. 

In order to render the Hamiltonian ([T]) more tractable, 
we make some general assumptions: Wc assume the elec- 
tron wave function in each dot takes the same form, 
thus allowing to be site-independent orbital energy levels 
£nia = £;<t, and Coulomb energies U n ji = Uji(j, I — d,b). 
We further consider the small dots to lie far apart from 
each other so that due to the large distance as well 
as the screening from electrons in the larger dots and 
background charges — which results in a small screen- 
ing length — the interdot Coulomb interaction for the 
localized states is small compared with their intradot 
Coulomb interaction. For simplicity, we shall consider 
only the Coulomb interaction between two neighboring 
b orbitals, i.e., Vbb(Rmn) = 5 m ^ n ±iU xbb where U xbb = 
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q—1 a, a' 
^ N N 
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where pda(q) = Hp=i d P+qa d piT an d the d-orbital en- 
ergy levels in momentum space read 
2tdCOs(27rp/N). We also define 



i(q) = E e-^"V tf (R n ). 



(13) 



Note that Vij{q) = since in the summation the term 
with n = is excluded; this term corresponds to the intra 
Coulomb interaction. For the lead-dot coupling t a d a in 
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Ht ©, after averaging out their site dependence, we 
have 



1 N 
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become 

+ H xU = Ubb n nb ^n nb l + U x bb 22 E 



N 



N 



n—1 o,a' 



N 



+ ^Wow(o) + ^ 2E"^^( )' 

n—1 crcr' 

where Uji = Uji+Vji(0)/2 andU xb b = Ua^exp [- Ri /L bb ] . 
The Hamiltonian comprising Eqs. ()2|9ll0ll7j) will be 
solved in Sec lIVI using a nonequilibrium Green function 
metho d 21 ' 22 . 



where k|| denotes the k component parallel to the plane 
of the QDA ring; rj denotes the additional degeneracy for 
the conduction electron state |k). We have absorbed the 
factor 1/y/N into t Q knb<r- 

To solve the Hamiltonian (|9I10I11I12[) would be ex- 
tremely difficult, since the Coulomb interaction between 
the d orbitals in momentum space is not diagonal. We 
now assume that the delocalized d states are in a para- 
magnetic phase; this phase is favored when tight-binding 
coupling is strong 2 ^, a condition taken throughout our 
numerical results. In this phase, the delocalized state in 
each quantum dot has equal probability of occupancy. 
Thus, on average (pda(q)) = Y, n (n n da)e lq ' t ' n equals zero 
unless q — 0. We then approximate pda{q) to be zero if 

Let us first look at the interdot Coulomb interaction 
between the delocalized states. In the paramagnetic 
phase, 
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5> dCT ((W(o). 
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We then can use Eq. (fl3|) to obtain 



Vdd(0) « U xd dj-j ^ 
n=l 



w ^ e [-NR lS in(im/N)/KL dd ] 



sin(7rn/JV) 



(16) 



for large N in a ID ring- like QDA. When Ldd — > oo, 
Vdd(0) ~ 2U x dd{l + ln2AT/7r), where 7 denotes Euler's 
constant. Hence, in the QDA of large size (N ~ 100), we 
find the total interdot Coulomb interaction to be compa- 
rable to the intradot Coulomb interaction (Udd) provided 
that we neglect the screening effect. Similar arguments 
hold for the interdot Coulomb interaction Vdb (0) between 
the localized and the delocalized states. For the interdot 
Coulomb interaction between localized orbitals, we shall 
keep only the nearest-neighbor term to examine its effect 
on the bistability. 

Combining the above results, the Coulomb interactions 



III. TUNNELLING CURRENT 

The general expression for the current is a function 
of the tunneling rate F and the local Green functions 
G r< , which are matrices spanned by the channels (or- 
bital and spin levels) 25 ' 26 . To start with, we assume there 
is no direct coupling between the localized and delocal- 
ized states, since their orbital wave functions have dif- 
ferent parities. Neglecting possible leaking inter-orbital 
tunneling current, T and G r< become block diagonal 
matrices with respect to the localized and delocalized or- 
bitals. The total current can then be divided into two 
parts, namely, I = lb + Id', h denotes the current tun- 
nelling through the I orbital (see Fig. [T] for the schematics 
of current tunnelling through a QDA.). 

We next assume the proportionality of the left and 
right tunneling rates for each current, so we can further 
rewrite the currents through the b and d states in terms 
of retarded Green functions only. They read^ 5 - - — 



de 



-w K 



lf L (e) - / fl ( e )]Im[tr(fGJ,)], (18) 



where T = T Le r m /(r Le + T m ); £ = b : d. f a (e) = f F {e- 
fi a ) and fi a are the Fermi distribution function and the 
chemical potential of the a-lead, respectively. W denotes 
the half band-width of the two metallic leads. Note that 
with localized states, we take the trace over the number 
of dots including spin, while with delocalized states, we 
take it over the band quantum number p(p = 1, • • • , N). 

For the localized orbitals, the current contributions 
from off-diagonal (inter-level) channels are of second or- 
der in tunneling rate T, since the off-diagonal Green 
functions G r n m for n ^ m (n, m labeling the chan- 
nels in the QDA) can be expressed in terms of diagonal 
Green functions multiplied by a self-energy S n>m , that 
is, G r n m ~ G^„E n , m G^ m ^. Moreover, neglecting the 
small real part of the self-energy, the imaginary part r, l m 
(Im£„ iTO ) reads 



TT Y2 tak.nbot a \fmba&{ £ — £ k) 
run ) 

2 k F R nut 



(19) 
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where t aknba = i Q! & ff e* k ir R ". T ba (e) = T Lb<J (e) + T Rba (e) 
with T aba (e)(= 2TrJ2 k \t a ba\ 2 S(e - e k )) being the tun- 
nelling rate from the \ba) state in the dots to the a-lead. 
Take the Fermi wave vector Uf ~ lOnm^ 1 for a nor- 
mal metallic lead with Fermi energy around 5 eV. Wc 
see that the factor s'm(kFR m n)/kFRmn is small when 
Rmn 3> lnm. Therefore, in the Coulomb blockade regime 
where the lead-dot coupling T is small, and with dots well 
separated, we can neglect those small off-diagonal terms 
in Eq. (fT5f and keep only the diagonal contributions. As- 
sume such to be the case for the localized states. We thus 
approximate the current I b to be 



h 



-Ne 



E 



^ J-w n rLbo-(e) + I 
x[f L (e)-f R (e)]lmGl a (e), (20) 



where G bcr (s) is the retarded 6-state Green function, 
which, by symmetry and degeneracy, is homogeneous for 
all dots. 

For delocalized states, from Eq. (fT4"|) we see that p is 
a good quantum number. Thus, the effect of coupling 
only leads to a shift in band energy, i.e. £dpa becomes 
£d P cr + RcSd + ird CT (e)/2. We then have 



Id 



N 



de T 



j-w 7T r Ldo .(e) + T RdtT (e) 

x[f L (£)-f R (e)}lmG;Ae), (21) 

where G r pcr {e) is the retarded d-state Green function. 

Because the band width of the metallic leads is com- 
monly much larger than the other energy scales in 
quantum dot systems, the tunnelling rates T a i (T (e) vary 
smoothly with energy. By assuming these rates to be 
energy- and spin-independent, we denote in a shorthand 
r a Zo-(e) = r Q i (I = d,b). To further facilitate analysis, 
we consider the system with a right-lead bias setting fi R 
lower than the local density of states of the QDA. In such 
case, we can safely set f R (e) = 0. We then find 



-Ne 



E 



ir 



deT Rb T Lb f L (e] 



ImGUe) 



NeT 



n ba ). 



(22) 



The tunnelling current I b and the &-orbital occupation 
number i^2 a {n br7 )) differ by a constant multiplicative 
factor eNT Rb /h. To study bistability in the tunnelling 
current, it suffices to compute the 6-orbital occupation 
number as a function of bias voltage. The tunnelling 
current Id and the average e?-orbital occupation num- 
ber (nda){= N^ 1 ^2p = i(n pcr )) can be related in a similar 
fashion, 



TVer 



Rd 



n do 



IV. GREEN FUNCTIONS 

To obtain the average 6-orbital and (i-orbital occu- 
pation numbers for tunnelling currents (|22|) and Id 
(|23|) . we used the noncquilibrium Green function method 
developed by Kuo and Chan g 21 ' 22 . This method keeps 
the self-energy from the lead-dot tunnelling coupling to 
lowest order, precisely corresponding to the scheme 2 of 
Ref. [29|. Similar to the Hartree-Fock approximatio n 30 ' 31 , 
we keep the self-energy term that describes the scatter- 
ing of a conduction electron by the electron in the dots, 
but neglect the spin-flip scattering processes responsible 
for the Kondo correlation s 30 ' 31 . The important differ- 
ence is that our method preserves the sum rule of the 
spectral density: — J d£lmG r /n — 1. Such approxima- 
tion scheme is generally considered valid in the Coulomb 
blockade regime, in which the lead-dot coupling is too 
weak for the Kondo effect to be observe d 29 ' 32 . Note that 
our lead-dot coupling in this Green function method is 
not limited to the weak-coupling limit, i.e.,. T k R T, 
which the standard Master equation methods generally 
assum e 33 ' 34 . The derivation of the Master equation in 
terms of the Keldysh Green function method has been 
reported 3 ^. For the reasons mentioned in the previous 
section, we also neglect contributions from off-diagonal 
elements. In a wide-band limit, the self-energy term gives 



£ 6 (e) = E 



ok 



*r b /2, 



E„(e) = V * iTdl 



(24) 
(25) 



where Ti = Tli + (I = d,b). The tunnelling rate Ti, 
being small, produces effectively a broadening width to 
the discrete energy levels for the quantum dots. 

It is interesting to note that the same Green function 
solution can be obtained by first solving the isolated QDA 
system and then replacing the positive infinitesimal imag- 
inary part 5 in frequency by a finite decay width T/2 in 
a phcnomcnological way^. 

We first derive the retarded Green functions. In the 
Zubarev notation^, a retarded Green function involving 
arbitrary fermionic operators A and B reads 



UA.B)) = -i lim 



dte i{z+iS)t ({A(t),B(0)}),(26) 



where (• • ■) denotes grand ensemble average and 6 a pos- 
itive infinitesimal. Using the Heisenberg equation of mo- 
tion, one can show 

{£ + i5)((A, B)) = ({A, B}) + (([A, H],B)). (27) 

This relation allows one to derive a flow of equations for 
the retarded Green functions. 

Since the \nba) states are homogeneous and degener- 
ate for all dots due to symmetry, we denote by G^(e) = 
(23) {(b a ,bl)) and G r bbrT (s) = ((n bs b a: bt)) the retarded b- 
state Green functions, with the site index n suppressed 



G 



(cr denotes the opposite spin of <r.); for retarded d- 
state Green functions, we specifically denote G prJ (e) = 
{(d pa ,$ p(T )), and G r ppry (e) = {{n^dp^d^)), since the 
\pda) band states are non-degenerate. We follow the pro- 
cedure of refs. Hi] and I22I First we obtain two Green func- 
tions for the \nba) state when all the \pda) band states 
(total 27V for TV dots) as well as the \ (n± l)ba) states are 
occupied: 



C^AT &ct( £ ) — 



(1 - Nbe)4Yl p Cp | 

fib — 4:U x bb — 2Udb fib — Ubb — 4:U x bb — 2Udb 

(28) 



G2N+1M ( e ) 



fib — Ubb — ^U x bb — 2Udb 



where Nbs = {nbs) is the single-particle occupation num- 
ber in the \ba) state. c p = (n pa .n pg .) is the correla- 
tion function between \pdcr) and \pda) band states, while 
c-b = {nbaTibs) is the correlation function between \ba) 
and \ba) states; fib = e — £ba + i^b/^- 

Next, the one-particle and two-particle Green func- 
tions for b orbitals are expressed as 



GUe) = (a' b + ti b + 



) 2 n^ 



G 



2N.bo 



n, 
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E E ft?" 
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n 



*66, i 
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fib — Ubb — n., 
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n 



db. 



(30) 



GlbAe) = {d'b + b'b + cbf Y[{d p + b p + c p ) 



5 3 

\ - \ - NbsPiPm 

z__l rn — 1 



n d6 , 



c ? Up Cp 

(31) 



The above operators and parameters are defined accord- 
ing to Refs. Hi] and The operators a p and b p are 
carried out explicitly in Eqs. (|30I31[) . Specifically, they 
are defined such that a p (f/g) = a p f/(g + 2Udb/N), 
bp(f/g) = b p f/(g + Udb/N), where 



a p = ((1 - n pa )(l - n p9 )) = 1 - N p(T - N pi ? + 
bp = (n pa (l - n pg )) + ((1 - n pa )n pd ) 
= N pa + N pg - 2c p . 



Similarly, the operators d' b and b' b are defined such that 

a' b (f/9) = abf/(g + 2U xbb ) and b' b (f/g) = b b f / (g + U xbb ). 
The parameters a pi b p , c p (a>b, 6&, Cb) are the probability 
weighting factors in the Green functions for the empty, 
singly-, and doubly-occupied configurations in the \pda) 
{\ba)) states, respectively. H x bb,i(= (i — ^)U x bb) denotes 
the interdot Coulomb repulsion felt by the \ba) state from 



its neighboring b orbitals in configuration i, while Tldb.m 
denotes the total Coulomb repulsion felt in configuration 
m, which depends on the number of electrons occupying 
the \pda) states — Physically, the Green functions pOI31[) 
compute the linear superposition of the probability am- 
plitudes of an electron propagating in the eigenstates of 
the isolated QDA system. To understand detailed deriva- 
tion of Green functions, readers are encouraged to see 
clear demonstration and discussion for the simple case 
of a nanostructure with two levels in Appendix. A of 
Ref. EI 
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The probability factors pt's are: p\ 
2a b b b , P3 = {b b + 2a b c b ), p4 = 2b b c b , pz = c b . In princi- 
ple, Eqs. (j3"0)) and (I3T1) allow one to compute the probabil- 
(29) ity factors p m in terms of the a p , b p and c p parameters for 
all 3 N terms via a simple recursion relation. In practice, 
when TV is large, the computation becomes impractical 
since it scales exponentially with TV. Fortunately, since 
we are considering a ring of quantum dots in which the 
\ba) states feel the same Coulomb energy Udb/N to all the 
\pda) states, we can reduce the summation of 3 N terms 
to 2TV+ 1 terms, because many terms will have the same 
denominator, and the computation scales linearly with 
TV. This observation immensely reduces the numerical 
complexity. The resulting Green functions read 



5 2N+1 



G r bAe) = J2 E 



PiPk 



i=l k=l 



fi b - n TbM - n' db k 

Nba 



fib — Ubb — H-xbb,i — II' 



db.k 



5 2N+1 



GlbAe)=J2 E 



NbaPiP'k 



i=l fe=l 



fib — Ubb — ^ X bbA — n^ fc k 



(32) 



(33) 



where II^ b k = (k—l)Udb/N and p' k is the sum of all those 
PmS that share the same Tl' db k , which can be obtained 
by using the following recursion relations: 

p' k (i) = a iP ' k (i - 1) + &ip' fe _i(i - 1) 

+CiP' k -2(. i - !); k < 2 «, 
P'S) = h iV' k -\( i -^) + c l p' k _ 2 (i-\)- 1 k = 2i 1 
p'S) = ap k _S - 1); k = % + 1, 

where p' k (i) denote the p' k coefficients at the i th iteration 
with Pi(0) = 1 and p' k (i) = for k < 1. The iteration 
process terminates at i = TV. 

With these retarded Green functions derived, we solve 
for the 6-orbital occupation numbers Nb a and Cb via 

Nha = j ^i ^^ ° b = J 2^t7 G "" t( ' £ - ) respec " 
tively, where the lesser Green function in the Coulomb 
blockade regime is related to the retarded Green function 
via 

TLbfL(e)+T R bf R (e) 



GfbM - 



-2i- 



-ImGL(e), (35) 
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We then iterate to self-consistency. Eq. (|35|) . while "ex- 
act" in the noninteracting model, also works well for the 
Anderson model in the Coulomb blockade regime. The 
derivation of lesser Green functions is similar to that of 
retarded ones. In [Appendix} a nanostructure system with 
two levels is given as an example to demonstrate Eqs. (|35|) 
and (|36|) by using the nonequilibrium equation of motion 
method^. To prove these two formulas for arbitrary 
levels can be done by the principle of induction, simi- 
lar to the way the retarded Green function is proved^. 
Rigorous proof will be given elsewhere. However, from 
the structures of retarded and lesser Green functions in 
the two-level system, as well as the fact that when all 
Coulomb interactions are turned off, the expressions re- 
duce to the exact solutions for noninteracting systems, it 
is not hard to see that Eqs. ()35I36|) hold true for arbitrary 
levels. 

At zero temperature, the average occupation numbers 
for the | bo) state read 



N, 



ha 



1 Lb 



5 2AT+1 



E E p*p'k 



(1 - N bs ) 



(37) 



x cot 



-N bs cot' 



i=l k=l 



r b /2 

x (V — £ba — Ubb — n K 66,i — n^ b k 
rV2 



r M 5 2JV+1 

-- — ^ — 2^ 1^ PiP k 

i=l fe=l 



7rr fi 



(38) 



x cot 



-1 



-6<T 



c^66 — n x 



66., 



IE 



db.k 



r b /2 



Besides some fixed physical parameters, it is worth not- 
ing that A^o- and c b are functions of the opposite-spin 
occupation number N bs and the applied bias V . 

Similar treatment can be used for the extended d-state 
occupation numbers N peT and c p . We have 



2JV-1 



i Ld 



E P' k (p){(l-N p *) 



(39) 



fc=i 



db cot 



1 fv Edpa n ddfe 



+66 cot 



-Cb cot 



-N„ 



T d /2 

_1 f V ~ £ dpcr ~ Udb - K'dd,k \ 
_1 ( V ~ £ dp<J -2Udb~ H'dd,k 
-1 (V — S dpa — ^dd,k+l 



-bb cot~ 



W2 

1 ( V - £ dpa- ~ Udb ~ Ii' d d,k+1 

{ r d /2 



+Cb cot 



V - e dV a - 2U db - n' 



cid,fc+l 



Td/2 



^LdNpa 



2N-1 



7rr d 

a& cot~ 



E f' fc (p) 



(40) 



k=l 

1 (V ~ S dp° ~ ^'dd,k+l 



-bb cot 



-Cb cot 



Td/2 

_1 ^ ^ - £ dpa - Udb - R'dd,k+1 \ 

Td/2 J 

-1 ( V ~ £ dp<? ~ 2U db - Tl' dd k+1 



r d /2 



where EE' 



dd,k 



(k — l)U~d d /N. p' k (p) is calculated ii 



ilar way to p' k except that the operator (a p + b p + c p ) is 
removed. The effect by the operator (a& + bb + c b ) is car- 
ried out explicitly when calculating G pa (e) and G ppa (e) . 
Note that a b and b b operators have different definitions 
from the a' b and b' b operators: a b (f / g) = a b f/(g + 2Ud b ) 

and bbif/g) = b b f/(g + U db ). 

To obtain the occupation numbers at nonzero temper- 
atures T, we first define a function -F(x, y) 



F(x,y) 



Im 



1 



V 



2TriksT 



, (41) 



where ^ denotes the digamma function. Thus, at 
nonzero temperatures, the average occupation numbers 
for the \ba) state read 



p 5 2N+1 

N *° = E pip* 
° i=i fc=i 



(1 - Nba) 

xF{ £brT + n X fe M + n' d6ifc ,r 6 /2) 



(42) 



+N bs F(s ba + u bb + n x66ii + k , T b /2) 



Cb = 



TLbNba 
7lT& 



5 2JV+1 



E E PiPk 
i=l fc=i 



x F(e 6CT + C/66 + n xb6)i + tt' dbJ; , r 6 /2) . (43) 
Similarly for A p(T and c p . 

V. RESULTS 



It is shown from Eqs. (f2"12j) and (|2"3"|) that the tunnelling 
current as a function of bias V can be related to the aver- 
age occupation numbers up to a constant multiplicative 
factor. It therefore suffices to investigate the average b- 
and d-orbital occupation numbers in order to understand 
the transport properties of the ring of coupled quantum 
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FIG. 2: (Color online) (a) The average b-orbital occupation 
number (nf, CT ) (b) The average d-orbital occupation number 
(rider) (= YLp=i ( n pcr) /iV) as a function of bias V (in units of T) 
as we sweep up (V t) and down (V 4-) the voltage bias. Differ- 
ent left-right tunnelling asymmetric factors 7?f sym = 1, 3, 9 
are taken, where i?f sym = TLd/^Rd- The curves at zero tem- 
perature T — show that the bistable effect enhances as 
the asymmetric factor 7?f sym increases.. We used N = 100 

dots, u = 7or, t d = 3or, E ba = e da = 3or, r 6 = 2r 

and Yd = 0.1F. The Coulomb energies are Ubt = U, 
Udt = Udd — 0.84(7, and U x bb = 0. All energy quantities 
are in units of a characteristic energy, T. 



dots. Wc shall consider two kinds of limiting cases: (i) 

tada <C t Q fccr and (ii) take <C t ac icr ■ 

The first kind can be realized in QDAs in which the 
smaller dot (which hosts the localized b orbital) is elon- 
gated, thereby allowing the orbital wave function to 
stretch out in the direction perpendicular to the plane 
of the ring (see Fig. [TJ. It then hybridizes with the elec- 
tronic states in the leads a lot stronger than the delo- 
calized orbital. This system will be used as an example 
in Sec IV AI Scc IV Gl to illustrate the effects on bistability 
produced by all the relevant physical parameters. 

The second kind can be realized in QDAs in which the 
small dots are embedded in a core-shell structure. It is 
natural for this case to assume that the tunnelling rate 
for the delocalized orbital is much larger. The tunnelling 
characteristics of the second kind and its difference from 



the first kind will be discussed in Scc IV HI 

Throughout the paper, all energy quantities are in 
units of a characteristic energy, T, which depends on sys- 
tems of interest. Furthermore, we take Tn, = = 
Tb/2, while considering various ratios of left and right 
tunnelling rates for the d-state, ^?f sym = ^Ld/^Rd, &s dis- 
cussed in Sec IV Al Varying the ratio T Lb /Tub produces 
similar effect. 

For the Coulomb interactions ([TT]) . we first set Ubb = 
U, which will also be used as a characteristic Coulomb 
energy. Wc then take interdot and intradot Coulomb in- 
teractions to be U X ji — Uji/2 and Udb = Udd = U/4, 
respectively. This will give Udb = Udd — U[l + (7 + 
ln27V/7r)/2]/4 without the screening effect. To clearly 
demonstrate the influences of the system parameters, 
particularly U and the QDA size N, over bistability in 
the QDA system, we will first focus on the ideal system, 
i.e., QDA without the inter 6-orbital Coulomb interaction 
and screening effect by setting U x bb — and Lji — > 00. 
These two effects will be considered in Sec IV Fl 

To understand bistability in the coupled dots systems, 
it is instructive to recognize its link with the population 
switching between the localized and delocalized states: 
Collectively, due to a large tight-binding parameter td, 
the d states form a finite band. The states in the lower 
part of the band get first populated as we increase the 
bias, see Fig. (Hd). Consequently, the energy level of 
the b states gets pushed up by the Coulomb repulsion to 
higher energy, therefore postponing the occurrence of its 
population, till those d states whose energies are below 
the effective energy level of the b states are all filled. 
Right as we continue increasing the bias when these d 
states are filled, a population switching between some 
higher-energy d states and the b states occurs. It results 
in a sudden jump in the tunnelling current as well as the 
average 6-orbital occupation number, and a sudden fall 
in the average d-orbital one. 

Reversely, as we decrease the bias from high-bias, we 
start to deplete the d states, which results in pushing the 
6-state energy level downward, till those d states whose 
energies are above the effective energy level of the b states 
arc all depleted, then a population switching occurs. 
Bistability shows up when these two events of switch- 
ing occur at different biases when we sweep upward and 
downward the bias voltage. Note, however, that the sig- 
nature of population switching does not guarantee bista- 
bility, as will be shown in Sec IV El We emphasize that 
the mechanism of population switching discussed here is 
to be distinguished in physical nature from that of the 
two-level quantum dots^ or the likes. 

While discussing in Scc IV Al and Scc IV El some of the 
results that were already shown by a simple mean-field 
approach^, we stress that our more complex method fo- 
cuses mainly on the temperature and size effects of the 
QDA system. Comparison between these two methods is 
made in Scc IV Bl 
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FIG. 3: (Color online) (a) The b-orbital occupation number 
(riba) (b) The average d-orbital occupation number (rida) as 
a function of bias V (in units of T) obtained by the mean-field 
theory and Kuo- Chang metho d 21,22 for a ring of coupled dots 
with N = 150, with -Rasym — 9. The other parameters are the 
same as in Fig. [2] 



Effect of asymmetry of left and right tunnelling 
rates 
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FIG. 4: (Color online) (a) The average &-orbital occupation 
number (riba) (b) The average d-orbital occupation number 
{nda) as a function of bias V (in units of T) for a ring of 
N coupled dots with N = 50, 100, 150, with Ri sym = 9. 
Udb = U dd ~ 0.75(7, 0.84(7, and 0.897 for N = 50, 100, 
and 150 respectively. As the size number N of dots increases, 
bistability enhances accordingly. Rest of parameters are same 
as in Fig. [2] 



From the above viewpoint of population switching, it 
is not surprising that the bistable effect is enhanced as 
we increase the d-orbital occupation numbers by increas- 
ing the degrees of freedom that form a band, like ex- 
tending from ID to 2D QDA or by tuning the left and 
right tunnelling rates. The latter effect on bistability is 
demonstrated in Fig. [3] in the limit t a d a <C t a \> ai where 
the asymmetric factor i?f sym = ^ Ldj^ Rd is varied. The 
figure shows that as i?f sym increases, so does the average 
d-orbital occupation number. As a result, the bistable 
effect enhances. 



B. Comparison with mean- field theory 

It is of interest to compare the solutions obtained 
from two different methods: the mean- field approachil 
and the Kuo-Chang metho d 21 ' 22 . It is well known that 



mean-field theory breaks automatically the particle-hole 
symmetry of the system Hamiltonian. Moreover, be- 
ing an effective single- particle approximation, it docs 
not treat two-particle correlation adequately, particularly 
when the correlation is strong. Recently, the noncquilib- 
rium mean-field theory has been criticized for possibly 
not being valid and giving false results, regarding bista- 
bility or hysteresis 2 ^. The present work, using a more 
careful treatment of finite-size quantum dot systems, val- 
idates in some respects the mean-field results that could 
otherwise raise doubt. 

Instead of treating two-particle correlation in a self- 
consistent way like the Kuo-Chang metho d 21 ' 22 , in mean- 
field theory, we decouple two-particle Green functions be- 
tween the localized b and delocalized d states, and those 
among the d states in the equations of motion into a prod- 
uct of a single-particle occupation number and a single- 
particle Green function. Namely, we make the following 
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approximations 



((n pa ,b a ,bl)) 



(n prT ,)((b a ,bt)), (44a) 
M((d pai dlj), (44b) 

(n pS )((V>4v))- (44c) 



Since we expect the Coulomb interactions between these 
states are suppressed by a factor of N (the total number 
of dots) , we expect the mean- field results to be sensible in 
the large N limit. Note, however, we leave {{nb S b a ,&£.)) 
intact for the localized b orbitals because its correlation 
is expected to be strong. Details of calculation can be 
found in Ref. [H 

Figure [3] compares the bistability results from these 
two methods. The results show that the mean-field ap- 
proach overestimates the bistable effect by 44%. The dis- 
crepancy for this effective single-particle theory is possi- 
bly attributed to artificial spontaneous symmetry break- 
ing and incorrect decoupling of two-particle correlations 
under the assumption that their correlations are small. 



C. Effect of QDA size 

A simple mean-field theory^ has demonstrated that 
devices made of ID and 2D QDAs can exhibit bistabil- 
ity in the tunnelling current as we tune the bias voltage. 
However, an important question remains as to how large 
a QDA must be in order to realize such bistability. This 
determines the density of memory bits made of such de- 
vices. The emergence of bistability in these systems is 
tightly related to the formation of a "good" band in one 
or two of the orbital states. But how do we judge if it 
forms a good band? In Fig. SI we show for large enough 
tight-binding parameter td and Coulomb energy U , the 
bistable regime (the bias interval where there exist two 
solutions for nonequilibrium states) starts to show up at 
N ~ 50. As we enlarge the size N of the ring of coupled 
dots (a ID QDA), the bistable effect enhances accord- 
ingly. This enhancement is caused both by a better band 
formation and by increasing Udb and Udd- 

The reason that a system of size N ~ 100 is enough to 
form a good (quasi) band can be understood by looking 
at the maximum energy level spacing for the d states, 
max({A p }) m Airtd/N and the effective Coulomb energy 

Udd/N felt among them. Figure [5] shows the spectral 
density for a ring of coupled dots with N = 100 as a 



function of energy at a bias before (see Fig. 5(a)) and 



after (sec Fig. 5(b) ) the population switching point (V 
71. ir 



Fig. 4(a) 



which can be seen from the dash-dotted curve in 
It is shown that the spectral density in the d 



states undergoes a tremendous redistribution before and 
after the switching point, while the peak in the spectral 
density of the b states shifts its position to the left. These 
sudden redistributions in the density of states reflect the 
sudden increase and fall of the average b- and <i-orbital 
occupation numbers, respectively. 
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FIG. 5: (Color online)The average 6-orbital spectral density 
Pba(s) = — ImGfKj and the average d-orbital spectral 
density pda (e) = — Im ~^2 p G prT (e) / (nN) as a function of fre- 
quency e (in units of V) at bias (a) V = 7ir (b) V = 71. 2T in 
the case of increasing bias, (a) and (b) correspond to the spec- 
tral density before and after the switching point V ~ 71. ir, 
respectively; see the dash-dotted curve in Fig. |4(a)| In plot- 
ting pdcris), we artificially broadened the spectrum by multi- 
plying T d by a factor of 20. We used R < L ym = 9, N = 100. 
The other parameters are the same as in Fig. [2] 



From the parameters used in Fig. [51 td — 30r and 
U = 70r, we get max({A p }) w 3.6r and U d d/N « 0.6r. 
On the other hand, the "effective" 6-state tunnelling rate 
Tb estimated from Fig. [5] gives roughly 6.4r and 7.4r 
before and after the switching point, both of which are 
a few times larger than the bare tunnelling rate Tb(= 
2r). With both d-state energy level spacings and effective 
Coulomb energy Udd/N being small, and a large effective 
&-state tunnelling rate, it is reasonable that, as far as the 
b state is concerned, it "sees" a good band, which triggers 
large population switching and consequently makes the 
bistable effect possible. 



D. Effect of tight-binding interaction td 

Here we investigate the dependence of bistability 
on the tight-binding interaction, td- The mean-field 
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FIG. 6: (Color online) (a) The average 6-orbital occupation 
number (n(, CT ) (b) The average d-orbital occupation number 
{iida) as a function of bias V (in units of F) for different tight- 
binding parameters td = 7r, 15r, 30r, AST for a QDA of size 
N = 100, with -Rasym = 9. Rest of parameters are the same 
as in Fig. [2] 
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FIG. 7: (Color online) (a) The average 6-orbital occupation 
number (n(, CT ) (b) The average d-orbital occupation number 
(rider ) as a function of bias V (in units of Y) for various 
Coulomb energy values U = 40r, 70r for a ring of coupled 

= 9. The other parameters 



dots with N = 100, with Rl 
are the same as in Fig. [2] 



approach^ concluded that the larger the td value is, the 
more effective bistability becomes. However, this is not 
the case when we consider a finite-size QDA. As was dis- 
cussed above, one necessary condition for bistability to 
take place is for the maximum energy level spacing among 
the d states to be sufficiently smaller than the effective 
6-state tunnelling rate. As shown in Fig. [5J for a ring 
of coupled dots with N — 100, increasing the td value 
will certainly bring d-state energy level spacings to be 
larger than the effective 6-state tunnelling rate (broaden- 
ing energy), thus destroying the (quasi-)continuous band 
structure and consequently bistability. For the four cases 
shown in Fig. [51 we estimate a suitable choice for td value 
to be around 15r(w E/ dd /4). 



E. Effect of Coulomb energy U 

A QDA of large enough size is not sufficient to realize 
bistability. Another condition is for the Coulomb energy 
U to be large. As is shown in Fig. [71 a small U does 
not exhibit bistability, even though there forms a good 
band in the d states. This is simply due to the fact 



that for a small U, the 6 state does not give the d states 
enough leverage to push its energy level sufficiently high 
to trigger bistability, even though population switching 
starts to show up. 



F. Effect of inter 6-orbital Coulomb interaction U x bb 
and of screening 

So far the inter 6-orbital Coulomb interaction and 
the screening effect are neglected. To illustrate these 
two effects, we take intradot Coulomb energies Ubb = 
U, Udb = 0.617, Udd = 0.5J7 and interdot Coulomb en- 
ergies U x db = Udb/4, U x dd = Udd/4:. We first study the 
role of inter 6-orbital Coulomb interaction over bistabil- 
ity. Next, we take into account the screening effect. 

Let us first consider the inter 6-orbital Coulomb inter- 
action and ignore the screening effect by setting Lji — > 
oo. In such case, we find Udb = 1.31Z7, Udd = 1.1J7, and 
U xb b = U xhh for QDA size N = 100. As shown in Fig. [5] 
for the 6- and d-orbital occupation numbers, inclusion 
of the Coulomb interaction between two neighboring 6 
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FIG. 8: (Color online) (a) The average 6-orbital occupation 
number (riba) (b) The average d-orbital occupation number 
(rider) as a function of bias V (in units of Y) for inter 6-orbital 
Coulomb energy values U x bb — and U/6 without screening 
(dash-dotted, solid curves), and U x bb = U/6 with screening 
(dashed curve). The screening lengths are Lbb = -Ri/2, Ldb = 
2Ri/3 and Ldd = Ri, where Ri is the nearest-neighbor dis- 
tance. We chose Ubb = U, Udb = 0.6C/, Udd = 0.5U and 
U x db = Udb/ '4, U x dd = Udd/ 4: (see the text for discussion). 
We also took the QDA size N = 100, U = 70r, t d = 20r and 
^asym = 9- The other parameters are the same as in Fig. [2] 



orbitals will decrease bistability. This suppression is re- 
flected on the sweep- down I — V curve, where the popu- 
lation switching point gradually shifts to the right as the 
interdot Coulomb repulsion strength U x bb increases. The 
switching point for the sweep-up I — V curve remains 
almost fixed. 

This phenomenon can be understood from the inter- 
play of probability factors pi's in Eq. ((57)) (In the consid- 
ered bias range, Cb is almost zero, and thus irrelevant): As 
we sweep upward the bias voltage, the 6 orbital stays un- 
charged, i.e., Nta ~ 0. This leads to ab »1, bb ~ 0, and 
Cb ~ 0, which leaves only the probability factor pi k 1, 
and the rest almost zero. This makes Eq. (|57|) and conse- 
quently the population switching point almost unaffected 
by U x bb- However, as we begin to sweep downward the 
bias voltage, the fact that Nb a ~ 1/3 at high bias makes 
the probability factors pi, pi and ps nonzero. The in- 
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FIG. 9: (Color online) (a) The average 6-orbital occupation 
number (nba) (b) The average d-orbital occupation number 
(rider) as a function of bias V (in units of Y) for temperatures 
k B T equal to 0, 3r, 6r for a QDA of size TV = 100, with 
Rasym = 9 and td = 15Y. The other parameters are the same 
as in Fig. [2] 



tcrplay of these three configurations (channels) leads to 
gradual decrease of Nba, instead of a sharp fall. Com- 
pared with the bistable range for U x bb = 0, it drops al- 
ready by 85% to 4.6T at U xbb = 11.67T(J7/6), showing 
that it is sensitive to the localized nature of the b or- 
bitals. 

To study the screening effect, we show how bistability 
is affected, through the suppression of interdot Coulomb 
interactions, by the values of screening length. From 
Eq. (|16[) . we find that the effective screened Coulomb 
energies Udb and Udd reduce only to 0.64J7 and 0.56U, 
respectively if we take the screening lengths Ldb = 2-Ri /3 
and Ldd = R\, where R\ is the nearest-neighbor distance. 
As discussed in Sec IV CI and Scc lVEl we expect that the 
reduction of Udb and Udd will suppress bistability. How- 
ever, as the impedimental inter 6-orbital Coulomb inter- 
action can also be screened to a much smaller value, we 
expect enhancement of bistability. For a shorter screen- 
ing length for the localized b orbital, namely Lbb = Ri/2, 
the screened inter 6-orbital Coulomb energy U x bb now re- 
duces to 0.l4U x bb- Thus, as shown in the dashed curves 
of Fig. |(5J), we still find a sizable bistable range 4.5T even 
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for small values of Udb and Udd- 



G. Finite temperature effect 

For practical applications, it is important to under- 
stand the temperature effect on bistability of this sys- 
tem. Figure [§] shows how bistability is destroyed when 
temperature increases. To explain the role of temper- 
ature, again we resort to the mechanism of population 
switching between the b and d orbital states. At zero 
temperature, when we increase the bias, the d states get 
first populated soon as their effective energy level meets 
the bias-dependent chemical potential, and then are able 
to push by repulsive Coulomb interaction the b state to 
higher energy. However, at finite temperatures larger 
than the tunnelling rates, the broadening effect due to 
finite tunnelling rates is taken over by the temperature 
which smears the characteristics of the spectral density. 
When the temperature equals the effective 6-state tun- 
nelling rate Tb, there is no priority for which of the d 
and b states to be populated. The population switching 
as well as bistability will therefore disappear. In Fig. |H1 
the rate Tb is equal to roughly 6.6r and 8.6r before and 
after the switching point (V = 68. 6r) at zero tempera- 
ture. Thus to ensure bistability, the temperature must 
be much smaller than Tb- 



H. The other limiting case: t a ba "C tada 

So far bistability is discussed by using the first case 
in which t a da t a ba- As for the opposite limiting case 
( t a ba -C t a da), we choose to show in Fig. [TD1 the tem- 
perature effect on bistability, to be compared with the 
corresponding results in the first case as shown in Fig. 1 

So far as bistability is concerned, there is no qualitative 
difference in these two opposite limiting cases. Quanti- 
tatively, however, the case in which t a b a *C t a d a shows 
a larger bistable effect. At zero temperature, its bistable 
range in Fig. [TU] measures 15.3r, larger than 11. Or in 
Fig. |U This quantitative advantage goes on to finite 
temperatures. However, at temperatures larger than 6r, 
bistability in both cases disappears. Generally, it shows 
that the second limiting case is more resistant to temper- 
ature, albeit within a characteristic temperature limit. 



VI. CONCLUSIONS 

In this paper, a noncquilibrium Green function 
metho d 21 ! 22 is applied to study the tunnelling current 
through a ring of coupled quantum dots (a ID QDA sys- 
tem). This method is able to treat adequately the elec- 
tronic correlation in nanostructures with multiple energy 
levels. Bistability in the tunnelling current is found if 
the number of dots is large enough and the temperature 
much smaller than the effective tunneling rates through 



0.3 
0.25 

0.2 
0.15 

0.1 
0.05 



. K B T=or,vl 
. K B T=or,vT 
. K B T=3r,vl 
. K B T=3r,vT 

K_T=6r,Vi 



60 

V 

(a) 




(b) 



FIG. 10: (Color online) (a) The average 6-orbital occupation 
number (ribo) (b) The average d-orbital occupation number 
{rider ) as a function of bias V (in units of Y) for temperatures 
k B T equal to 0, 3r, 6r for a QDA of size N = 100. In 
contrast with all the other figures, here the opposite limiting 
situation is considered, namely, we take Tb = 0.2r, Fd — 2F 
The other parameters are the same as in Fig. [9] 



the localized states. The effects due to size scaling, 
Coulomb energies, screening effect, tight-binding inter- 
action for the delocalized states on bistability were also 
investigated. Our results show that in order for bistabil- 
ity to occur, we need to have a) a large enough Coulomb 
energy U, b) small intcrdot Coulomb interaction (U x bb) 
between the localized states, and c) the tight-binding in- 
teraction between the delocalized states (td) to be around 
£7/4 for best performance. Furthermore, we showed the 
mixed role of the screening effect: on the one hand, it 
decreases bistability through reduction of the effective 
Coulomb interactions Udb and Udd] on the other hand, 
it enhances bistability through reduction of the inter b- 
orbital Coulomb interaction U x bb- Most importantly, we 
show that bistability in the tunnelling current occurs 
when ./V > 50 and becomes sizable around N ~ 100 in 
a QDA system. This size estimation provides an impor- 
tant information for fabrication of QDAs aimed at mem- 
ory devices. Based on our analysis for a QDA of size 
TV = 100, in order to use the ID QDA as a memory de- 
vice at room temperature, the characteristic energy (r) 
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used must be around 10 meV and a Coulomb charging 
energy around 700 meV, which can only be achieved in 
very small quantum dots (with size around lnm). The 
requirement should become less stringent for ID QDAs 
of larger size, and for 2D QDAs as demonstrated in the 
previous study^ based on mean-field theory. 
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Appendix. LESSER GREEN FUNCTIONS FOR 
A TWO-LEVEL QUANTUM DOT 

In this appendix, we shall derive in detail the lesser 
Green functions ()35I36|) for the case of a two-level nanos- 
tructure using the same approximation mentioned in 
ScclIVI Consider two levels denoted by I and j (j ^ I). 
Using the nonequilibrium equation of motion method^, 
the equation of motion for the lesser Green function 
Gf a {e) reads 

W G< CT ( £ ) = E< (e)G^(e) + U t Gf/e) 

+U e . J [G<^(e) + Gf^Js)], (A.l) 

where \n = e — Ee a + iTg/2. The lesser self-energy 
is given by £<(e) = i J2 a =L,R ^odcfa^)- G£ CT ( £ ) 
denotes the advanced Green function. Gf e (e) = 

((fi£,S*, t 7,4, t r)) < > G f,j,a( £ ) = (( n j,* d Z,<r> d l,<r)) < and 

Gfj (e) = ((rij t( rdg ta , d\ <J )) < are two-particle Green 
functions coupled to Gf a via the intralevel and interlevel 
Coulomb interactions. 

To solve Eq. (|A.lj) . we derive the equations of motion 
for the two-particle Green functions. They satisfy for 

[m-Ut)Gf,{e) = Vt,(e)C%Ae) + V t j\Gf, i J S ) 



+Gf,iJe)], 



(A.2) 



( W - U u )G<{e) = S< ( £ )G? , Je) + U ( G<(e) 



+U ( . ] G<(e), 



(A.3) 



On - U u )G<^{e) = E< (e)G^- a (e) + U e G< Aj Je) 
+U e ,jG< jtj (e). (A.4) 

Two-particle lesser Green functions arc now 

coupled to the following three-particle Green 

functions: Gf g j -(s) = ((nt^rij^de^, 4, < r)) < ' 

G tej,<r( e ) = (( n e,v n j,<r d t,<r,dlv)) < , and Gf dj (s) = 



((nj^rijicrde^, d\ cr )) < - These three-particle Green func- 
tions will again couple via their equations of motion with 
the four-particle Green functions, where the hierarchy 
terminates. Three-particle Green functions are given by 

(»l-U e -U itj )G< Aj ^(e) = Ef a (e)Gl e ^(e) 

+U e ,jGf(e), (A.5) 

(pi-U e -U £j )G< e ^(e) = £<( £ )G? /iii<7 ( £ ) 

+U e ,jGf(e), (A.6) 
( W -2U^G<..( £ ) = X<(e)Gi^(e) 

+U l Gf{e), (A.7) 

where Gf(e) = (X n i,a n 3,a n j-<J d i.<7-> d \ cr )) < is the four- 
particle Green's function whose equation of motion reads 

{ lH -U t -2U^)G<{e) = ££(e)G2(e). (A.8) 

To solve the closed set of lesser Green functions requires 
solutions of the advanced Green functions G°'s. Their 
complex conjugates, namely the retarded Green func- 
tions, can be found in Appendix. A of Ref. With 
G%{e) known, we immediately see 



G 4 <( £ ) = E<(e) 



Nt 9 Ci 



\m-u t -2U id \ i 



(A.9) 



Inserting the above equation and solutions of G"^ ^ 

G t,i,j A e ) and G ljj& int0 Ec l s - <|A.5IA.6IA.7j) yields 
three-particle Green functions as 



G f,e,j,cr( £ ) 



+- 



',z)Nt,a 

c 



\m-u e -u itJ f 



U £ -2Ut tj \ 2 



(A.10) 



^•j(e) = E^(e)cj 



1 - N f 



■\lH-2U id \ 2 



\ f H - Ue - 2U Lj \ 2 



(All) 



and 



2c, 



\fj,i -U t - U e J< 



\Hi-Ut-2UijY 



(A.12) 



where bj, Cj as well as aj are defined in the main text. 
Inserting Eq. (|A.12[) and the solution of G" f ( £ ) into 
Eq. |A~2|) yields 



\Hl-U t \ 2 \iH-Ut-U ttj \ 2 



(A.13) 
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while inserting Eqs. (|A.11[) , (|A.12[) and solutions of where III = 0, II2 = Ut,j, and II2 = lUi.y, p\ 



G a i]S (s), G1 ] a {e) into Eqs. fO]) and (|A~4"]) , we obtain 

1 - N,s 



Xe){l 



\ f M -U t - C/^f 



2Cn 



\»i-Ui,j\ 2 



\H-2U td \ 2 



\lH-Ui-2Utj\ s 



}.(A.14) 



Combining Eqs. (|A.13[) . (|A.14|) . (|A.1[) and the solution 
of Gg a (e), we obtain finally 



m— 1 

+ 



(1 - jV) 
lw-n m | 2 



m - Ut- n m | 2 



(A.15) 



p 2 = bj, and p 3 = cj. 

From the expressions for all these lesser Green func- 
tions, we see that they all take the form of Eqs. (|35l36p . 
namely, 



= -2^MlmG'' 

= _ 2 , r ^-fc( £ )+ r ^/«( £ ) ImG r (A . 16) 
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